Vibrational Andreev bound states in magnetic molecules 



(N 

o 

(N 



i 

G 

O 

o 



> 
in 

On 

in 
O 
(N 



X 



Denis Golez, 1 Janez Bonca, 1,2 and Rok Zitko 1 - 2 

'jozef Stefan Institute, Jamova 39, SI- 1000 Ljubljana, Slovenia 
2 FacuIty of Mathematics and Physics, University of Ljubljana, Jadranska 19, SI-1000 Ljubljana, Slovenia 

We predict the existence of vibrational Andreev bound states in deformable magnetic molecules on super- 
conducting surfaces. We discuss the Anderson impurity model with electron-phonon coupling to a realistic 
anharmonic vibrational mode that modulates the tunneling barrier and show that the vibronic features are spec- 
troscopically most visibile near the transition point between the Kondo-screened singlet and the unscreened 
doublet ground state. We find competing tendencies between phonon hardening due to anharmonicity and soft- 
ening due to coupling to electrons, contrary to the Anderson-Holstein model and other models with harmonic 
local phonon mode where the vibrational mode is always softened. In addition, we find that the singlet and 
doublet many-body states may experience very different effective phonon potentials. 

PACS numbers: 72.15.Qm, 73.20.Hb, 73.40.Gk, 75.75.-c, 74.55,+v 



Molecules that conduct electrical current when em- 

bedded in a junction between two metal electrodes fl-can 
become the active element in a circuit, such as a rectifier 16j|7|] 
or a memory element . Alternatively, molecules deposited 
on a metal substrate can be probed by a scanning tunneling 
microscope to study their diffusion fl, conformation changes 
JToL [Till , dissociation 1 121. and chemical reactions (13, 14]. 
Strong coupling between electronic and vibrational degrees 
of freedom plays a critical role for the molecule's functional 



properties H151 I16H . The electron-phonon (e-ph) couplin 



renormalizes the electron-electron (e-e) interaction [|17l 
and, if large enough, leads to an effective attractive interaction 

112, 
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2011 . Vibrational modes are detected in the differential 
conductance spectra as spectral features at characteristic fre- 
quencies I21M24I1 . that serve as "molecular fingerprints" 1125 - 
2711 . In magnetic molecules on normal metal surfaces, the low- 
temperature spectra exhibit zero-bias anomalies due to the 
Kondo screening of the local moment 1 28, 29j], while magnetic 
molecules adsorbed on superconductors exhibit sharp spec- 
tral peaks inside the gap (Andreev bound states, ABS) due 
to the c omp etition between the Kondo effect and the electron 
pairing jl30l - l35ll . as has been recently experimentally demon- 
strated 13611 . Since molecules are deformable, the vibrational 
modes need to be taken into account for a comprehensive de- 
scription of all features that may occur in the subgap part of 
the spectrum. 

In this work we study a realistic model of a deformable 
magnetic molecule in contact with a superconductor. The 
vibrational mode is described using the anharmonic Morse 
potential 11371 13811 and the displacement exponentially modu- 
lates the tunneling barrier [39]. The anharmonicity is required 
to remove unphysical infinite-displacement solution found in 
the harmonic approximation, while the exponential modula- 
tion removes the fluctuating-sign problem of the lowest-order 
linear coupling; both choices are also closer to reality. We 
focus on the case of a weakly bound adsorbate with external 
(center-of-mass) vibrational mode whose energy is compara- 
ble with the superconducting gap, so that vibronic features 
occur inside the gap. In particular, we show that the molecule 
spectral function features vibrational side-peaks in addition to 



the main ABS peaks. The side-peaks are visible for generic 
model parameters even for moderate, experimentally relevant 
e-ph coupling strength. Because the peak width is due only to 
thermal broadening and experimental noise, this setup permits 
very precise determination of the phonon frequency renormal- 
ization due to the e-ph coupling and the e-e interactions. 
We describe the system with the impurity model H = 

-ffband + H mo \ + H osc + H coup . Here i?band = 

Ek* e kc{* c k,v + A Efe, <J ( c L ct -fe,- ( T + h - c -)> where C k,° 

are the conduction-band electron operators with momentum 
k, spin a and energy e^, and A is the superconducting gap. 
Hmo\ = e(n-|- + n\) + Un-fTii is the molecule Hamiltonian, 
where n c = S a d a , e is the on-site energy, and U the e-e re- 
pulsion. The quantity 5 = e + U/2 measures the deviation 
of the system from the particle-hole (p-h) symmetric point. 
The displaced molecule feels a realistic Morse potential of the 
form 



VWsc = D e [1 - exp (-bx)} 



(1) 



where D e is the well depth and b controls its width. The oscil- 
lator Hamiltonian is thus H osc = p 2 /2m + VMorsc- We define 
the harmonic frequency as uq = by/2D e / m, and the displace- 
ment and momentum operators as x = (a + a^)yl/2mwo 
and p = i(a' — a) where d and d^ are the phonon 
ladder operators. The shape of the potential is then fully de- 
scribed by two parameters, D e and ojq. The harmonic poten- 
tial is recovered in the limit D e —> oo, keeping u>o constant. 
In the coupling part H coup = V(x) J2k.a( c ka d °- + n - c -)> the 
tunneling term is exponentially modulated by the molecular 
displacement: 

V(x) = V exp [-g(a + d f )] , 

where g > is the e-ph coupling constant. The hybridization 
strength at zero displacement is characterized by r = ttpVq, 
where p is the density of states in the band in the normal state. 
In all numerical calculations presented in this work, we set 
U = 1, A = 0.04, and cu = 0.01. 

In the absence of e-ph coupling (g = 0) and close to the 
particle-hole symmetric point (S ~ 0), the ground state of the 



2 



system depends on the relative values of the Kondo tempera- 
ture T K and the BCS gap A JM^l. In the limit T K > A, 
the local moment is screened by non-paired conducting elec- 
trons. The ground state is then a spin singlet (S) many-body 
Kondo state (3(itl . In the opposite limit of Tk <C A, the forma- 
tion of the Cooper pairs leaves no low-energy electrons avail- 
able to screen the impurity spin. The ground state is then a 
spin doublet (D) unscreened many-body state. For Tk ~ A, 
there is a level crossing between the two different ground 
states. The subgap part of the spectral function for g = 
has generically two peaks, located symmetrically with respect 
to the Fermi level (here fixed at oj = 0). These features, of- 
ten referred to as the Andreev bound states or Shiba states, 
correspond to the transitions between the S and D many-body 
states. These states are "bound" in the sense that they cor- 
respond to excitations below the continuum of quasiparticle 
states and that the corresponding spectral weight is localized 
around the impurity site. The transitions can occur either by 
injecting an electron (cj > peak) or by removing it (uj < 
peak). Away from the p-h symmetric point, the weights of 
these two peaks are different B35I1 . Their spectral weights go 
to zero as the S-D energy difference tends toward A in both 
small-r (T K < A) and large-r (T K > A) hybridization lim- 
its and is maximal near the S-D level crossing B4I1 . which is 
signaled by the crossing of the subgap peaks in the spectrum 

We now consider the effect of the e-ph coupling (g > 0) 
on the subgap states. For small coupling the hybridization 
is renormalized as Tq — >• T — T(g, (x)), where (x) is the 
expectation value of the displacement that is negative (the 
molecule approaches the surface) and linear in g [see Fig. 3(g)- 
(i)], thus the enhancement of the hybridization is approxi- 
mately quadratic. For large coupling this mean-field approx- 
imation is no longer accurate, but the trend toward stronger 
effective V remains. In addition, if the phonon frequency is 
smaller that the BCS gap, luq < A, entirely new features 
(vibronic ABS sidepeaks) appear in the subgap part of the 
molecule spectral function, as we show in the following. 

The full prob lem is solved using the numerical renormaliza- 
tion group 1 40l 41 ] with extensions for superconducting bands 
This technique treats all interactions (e-ph coupling, 
e-e repulsion, BCS pairing) on equal footing. The calculations 
have been performed for the discretization parameter A = 4 
and are fully converged with respect to the phonon cutoff and 
the number of states kept in the truncation. 

We first consider the evolution of the subgap states as a 
function of the hybridization strength To, keeping the e-ph 
coupling g constant, see Fig. da). This is motivated by the ex- 
perimental realization of the singlet-doublet crossing, which 
is tuned by small changes in the molecule-substrate coupling 
J36ll . The calculation has been performed for generic parame- 
ter values, in particular away from the p-h symmetric point. 
The level crossing between the doublet and singlet ground 
states DO and SO occurs at Tq = Tq c - All other excited sub- 
gap states, one series of D states (Dl, D2, . . . ) and another 
of S states (SI, S2, . . . ), are phonon induced. For a harmonic 
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FIG. 1. (Color online) (a) Energy of the subgap singlet (S) and dou- 
blet (D) states as a function of the hybridization strength To. (b) 
Subgap part of the spectral function. Subfigures show cross-sections 
at particular values of To, indicated by dashed lines, (c) Visibility of 
the first and second vibronic side-peaks. Parameters are g = 0.05, 
D e = 0.5, and 5 = 0.1. 



potential and very weak e-ph coupling, these vibronic states 
would form ladders with equidistant spacing ljq. For the an- 
harmonic Morse potential, however, the exact solution for the 
vibrational energy levels includes a quadratic correction term 



I 37H and, furthermore, the effective e nerg y spacing is renor- 
malized by the e-ph coupling see also Fig. |3jd)- 



(f). 

The subgap part of the spectral function, Fig. [Tib), demon- 
strates that in addition to the main spectral peak (S0-D0 transi- 
tions), other phonon-induced side-peaks (S0-D1, D0-S1, etc. 
transitions) are also present and have sufficient spectral weight 
to be observable. The subfigures Fig. [JJbl-b3) show spec- 
tral curves at three characteristic parameter regimes; first side 
peaks are clearly resolved in all cases. The ratio between the 
spectral weight of the side peaks and that of the main ABS 
peaks (i.e., the visibility), Fig.[TJc), is maximal near the level 
crossing. We also find that the visibility is largest at the p-h 
symmetric point (S = 0) and is reduced somewhat away from 
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FIG. 2. (Color online) (a) Energy of the subgap singlet (S) and dou- 
blet (D) states as a function of the electron-phonon coupling g. (b) 
Subgap part of the spectral function. Parameters are To = 0.1, 
D e = 0.1, and 5 = 0. 



it (results not shown). 

Alternatively, the S-D transition can be induced by increas- 
ing the e-ph coupling g at fixed To, see Fig. |2] The vibra- 
tional excited states exhibit an unexpected feature: with in- 
creasing e-ph coupling, the effective phonon frequency in- 
creases (the phonon mode hardens). The renormalization of 
the phonon frequency has been noted in p revious studies of 
impurity models with vibration modes ]39j |42|444J| , where the 
phonon mode softens. This is the case both for the coupling 
to charge (Anderson-Holstein model) and for the coupling to 
the center-of-mass modes. The phonon hardening is a char- 
acteristic feature of the anharmonic potential and results from 
the displacement of the oscillator to the part of the potential 
with higher second derivative, while the softening results from 
the electron tunneling fluctuations l45ll and occurs generic ally. 
Both tendencies are present in our model and the dominant ef- 
fect depends on the value of the parameter D e . 

In order to compare the degree of the anharmonicity, the 
bare (non-renormalized) potential profiles are shown in the 
top panel of Fig. [3] In Fig. 0a)-(c) the energy dependence 
of the Andreev states shows the S-D transition as a func- 
tion of the e-ph coupling g for different values of D e . The 
renormalized phonon frequencies, Fig. |3jd)-(f), indicate that 
for D e = 0.1, the potential is strongly anharmonic and the 
phonon mode hardens, while for D e = 5, in the harmonic 
limit, it softens. We also observe remarkably large differ- 
ences in the D and S sectors, which are the most pronounced 
for intermediate anharmonicity where, as a function of g, the 
phonon mode softens in the D sector, while it hardens in the 
S sector [see Fig. HJe)]. The deformation of the molecules 
as a function of g is shown in Fig. |2g)-(i) [notice that the 
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FIG. 3. (Color online) Top panel: bare (non-renormalized) Morse po- 
tential for different parameter of the anharmonicity. Bottom panel: 
(a,b,c) Energy of the subgap singlet (S) and doublet (D) states ver- 
sus electron-phonon coupling g for different parameter of the anhar- 
monicity. Parameters: ro = 0.1, 5 = 0. (d,e,f) Effective phonon 
frequencies of D and S states. (g,h,i) Displacement and fluctuations 
(insets) of the displacement. (j,k,l) Visibility of first vibronic side- 
peak. 



horizontal ranges are different]. In the harmonic limit the e- 
ph coupling leads to unphysically strong deformation already 
at small values of g, while for anharmonic potential, the dis- 
placement of the molecule is constrained by the repulsive part 

of the Morse potential. The insets represent the fluctuations of 

l Ii 

the displacement, Sx = {(x 2 ) — (a;) 2 ) , which in the har- 
monic limit strongly increase at the transition and then rapidly 
drop at higher g, while for strongly anharmonic potential they 
decrease monotonously and there is no enhancement at the 
transition. Independent of the anharmonicity parameter D e , 
the visibility grows with increasing e-ph coupling for g < g c , 
see Fig. [3]j)-(l). For larger e-ph coupling the visibility starts 
to decrease, since the energy difference between the ground 
state and the side-peak states approaches the BCS gap 134 1. 
The maximum of the visibility thus genetically occurs near 
the S-D transition. The visibility is larger for high D e values, 
i.e., in the harmonic limit. 

The oscillator distribution functions p(x) and the effective 
potentials V e s(x), calculated using the reduced phonon den- 
sity matrix |46|], are shown in Fig. |4] The minimum of the 
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effective potential is significantly shifted due to the e-ph cou- 
pling to its new equilibrium position and we find that the 
effective potential is not the same for S and D states. This 
implies that the oscillator parameters are renormalized differ- 
ently in the two spin sectors due to different electron tunneling 
rates in the unscreened D and Kondo screened S many-body 
states. In the weak coupling (small-g) regime the effect of the 
phonons on the subgap part of the spectral function can be un- 
derstood using the Frank-Condon principle |47[, which states 
that the transition between the vibrational states is more likely 
to happen if the vibrational wave functions overlap more sig- 
nificantly. The overlap of the effective vibronic part of the 
wave functions, 4>{x) = \J p(x), is proportional to the transi- 
tion probability between the states. For the lowest-lying sin- 
glet SO and doublet DO states, and the excited singlet 51 and 
doublet Dl states, the overlaps are represented in Fig. |4}d). 
Within the Franck-Condon approximation the overlap ratio 
\(S0\DI)\ 2 /\(S0\D0)\ 2 is proportional to the visibility and 
provides good agreement with the actual values, Fig. |3jj), 
for small electron-phonon couplings g. However, since the 
Franck-Condon principle is based on the Born-Oppenheimer 
approximation, polaronic effects are neglected and, therefore, 
it cannot explain the features in the spectral function for e-ph 
interaction close to or beyond the level crossing. The sud- 
den increase in the overlap between the SO and Dl states for 
g « 0.115, coincides with the transition of the Dl state into 
the continuum, where the excited Dl state strongly mixes with 
the direct product states consisting of the SO state and an addi- 
tional continuum low-energy quasiparticle, thus its vibrational 
properties are essentially those of the SO state, which explains 
the perfect overlap. 

We have shown that due to the deformability of the 
molecule additional vibronic states occur inside the BCS gap, 
which are clearly visible in the subgap part of the spectral 
function. We find that generically the intensity of the side 
peak reaches a maximum close to the ground-state level cross- 
ing, but that the peak intensity depends on the system param- 
eters: it increases near the particle-hole symmetric point, for 
large D e , and for large e-ph coupling g. The increase in the 
visibility for weak e-ph coupling can be described using the 
Franck-Condon principle. 

To test the predictions of this work, we propose as candi- 
date systems planar macrocyclic molecules which form coor- 
dination complexes with weak metal-ligand bonds, so that the 
magnetic ions support low frequency "rattling" modes, while 
appropriate surfaces are elemental BCS superconductors. 

The authors acknowledge discussions with Jernej Mravlje, 
Tomaz Rejec, Anton Ramsak and Peter Prelovsek, and the 
support from the Slovenian Research Agency (ARRS) under 
Program PI -0044. 
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FIG. 4. (Color online) (a) Effective (renormalized) potential of S 
and D subgap states and bare (non-renormalized) Morse potential. 
(b,c) Oscillator distribution p(x) for ABS in the effective potential 
for S and D states. The probability is not normalized, (d) Overlap of 
the phonon part of the wavefunction between different subgap states. 
Parameters are T = 0.1, D e = 0.1, 8 = 0. 
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